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Abstract 



o 



We employ a fast and accurate algorithm to treat dielectric interfaces within molecular dynamics simulations and demonstrate the 
importance of dielectric boundary forces (DBFs) in two systems of interests in soft-condensed matter science. We investigate a salt 
solution confined to a slit pore, and a model of a DNA fragment translocating thorugh a narrow pore. 
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1. Introduction 

Coarse-grained models are widely employed in computer 
simulations of soft-condensed matter systems because of the 
significant computational speed-up granted by the reduction in 
the number of degrees of freedom. A common coarse-graining 
approach is to employ an implicit solvent approach, thereby 
removing completely all solvent molecules. In the simplest 
scheme the polarizability of an aqueous solvent is taken into 
account by setting the dielectric constant to 80, thus reducing 
the electrostatic interaction by the same factor. This approach, 
however, fails in presence of interfaces between the solvent and 
materials of significantly different polarizability. Various meth- 
ods to solve the Poisson equation with inhomogenous dielectric 
permittivity are available, and have been employed in different 
contexts, such as in Refs. [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1 1]. One of 
the possibilities consists in solving the boundary integral equa- 
tions using boundary element methods. This approach is partic- 
ularly useful when the boundary integral problem is formulated 
in terms of an system of implicit linear equations which express 
the induced surface polarization charges in terms of the elec- 
tric field[ll, 12, 13, 14, 15, 16]. With the ICC* algorithm[17] 
(Induced Charge Calculation with fast Coulomb Solvers) we 
enhanced the capabilities of different, widely employed elec- 
trostatic solvers with an iterative scheme to solve this set of 
equations. Noticeably, the procedure automatically yields so- 
lutions which obey the boundary conditions of the underlying 
Coulomb solvers, which is important to reduce the necessary 
system size. 

In this article we review the algorithm and present two appli- 
cations in the context of soft-condensed matter. The first one 
is the test case scenario of an electrolyte confined between two 
walls. In the second application we calculate the free energy 
barrier which a model DNA fragment has to overcome in order 
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to be transported through a synthetic nanopore, for two different 
ionic strength conditions. 

2. The ICC* Algorithm 

The goal of the ICC* algorithm is to solve the Poisson equa- 
tion for an inhomogenous dielectric. The Poisson equation in 
cgs units reads as: 



V • (eV<D) = -4np ext , 



(1) 



Let us suppose for simplicity that there is only one interface 
between the two regions 1 and 2 of different dielectric permit- 
tivities ei and E2- In order to express the boundary conditions 
in terms of induced surface charges, we integrate eq. 1 over a 
pillbox small enough that the inhomogenity of the electric field 
through its caps is negligible. Then the induced charge in the 
pillbox can be expressed as a surface integral 
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(2) 



where n is the surface normal, conventionally pointing from 
region 1 into region 2, and Ei/a is the electric field in medium 
1 or 2 in close proximity of the interface. The last equation is 
valid in the limit of a surface element of infinitesimal area A. 

The field E1/2 at a given position on the surface can be written 
as a superposition of the field generated by the charged surface 
element in the very same location and of the field generated by 
other external and induced charges in the system. We denote the 
second contribution by E and the charge density on the surface 
element by <x = qi„d/A. 



E1/2 = E + 2n/sicm. 



(3) 



By combining Eqns.(2) and (3) one obtains an expression for 
the induced surface charge: 
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This boundary integral formulation can be implemented in a 
numerical boundary element scheme by assigning charge ele- 
ments on a grid which discretizes the interface. Eq.(4) then 
is a set of coupled equations, which can be solved with dif- 
ferent approaches. In the ICC* implementation we decided to 
use a successive over-relaxation (SOR) scheme to solve them 
iteratively. When all boundary elements are approximated as 
point charges it is possible to determine E with widely used 
fast Coulomb solvers that take into account the desired periodic 
boundary conditions of the system, and simultaneously deter- 
mine the surface charge density to arbitrary precision. In every 
step of the iterative scheme the new guess for the charge at each 
surface element q' nev/ is determined from the previous value with 
the following equation: 



(5) 



where <x, is obtained from Eq. (4). Here A is a free parame- 
ter in the range between and 2. With a value of A = 0.9 
we obtained a satisfactory speed of convergence and no stabil- 
ity issues in all performed simulations. Although more refined 
approximations than the use of point charges are available [10], 
our choice allowed us to readily implement the ICC* algorithm 
in the ESPResSo molecular dynamics package. In ESPResSo 
it is currently possible to use ICC* with the following electro- 
static solvers: P3M[18, 19, 20], ELC[21, 22], MMM2D[23], 
and MMM1D[24, 25]. 

Finally it should be noted that in a typical coarse-grained 
simulation, the small change of particle positions in one sim- 
ulation timestep leads only to a small change in the boundary 
element charge, so that each ICC* update usually needs only 
1-3 steps of the iterative algorithm. 

3. Electrolyte in a Slit Geometry 

As a first example of the use of ICC* we report here the re- 
sults of a coarse-grained Langevin dynamics simulation of a 
salt solution confined to a slit pore. We employ the restricted 
primitive model for electrolytes, namely representing ions as 
repulsive charged spheres embedded in a dielectric continuum 
with s = 80, and use the ICC* algorithm to investigate the role 
of dielectric contrast between the solution and the wall material. 
The temperature T was set to 300 K, and the Bjerrum length to 
0.71 nm, corresponding to that of water at room temperature. 
An excluded volume interaction between ions was set up by 
means of a Weeks-Chandler- Anderson potential 



Uu(r) = 



4e, 



otherwise, 



(6) 



with interaction strength eu equal to IkgT {kg being Boltz- 
mann's constant) and a cutoff radius r c - 2^ 6 cr, where <x - 
0.284 nm. The slit pore has a dimension of 7 x 7 x 3.5 nm, 
with periodic boundary conditions applied along the two longer 
edges, and has been filled by 50 positive and 50 negative ions, 
reaching an approximate concentration of 200 mmol/1. In order 
to investigate the effect of the dielectric mismatch, the dielectric 
permittivity on the left wall of the pore has been set to 6400, 




1.5 2 
Position / nm 

Figure 1 : Density profile of ions in a slit geometry. The dielectric constant 
of the left, middle, and right regions are 6400, 80 and 80, respectively. The 
ICMMM2D (solid line) and ICC* results (squares) are reported. 



while no dielectric contrast has been set for the solvent/right 
wall interface. In Fig.l we present a comparison between the 
density profile of ions in the slit, obtained with ICC* in com- 
bination with MMM2D as an electrostatic solver, and with 
ICMMM2D. The latter method is an extension of the MMM2D 
algorithm, that takes into account the presence of dielectric mis- 
match for flat boundaries by means of image charges. As it is 
apparent from the results in Fig.l, the use of a discretization 
grid of 10 x lOx elements for ICC* is enough to obtain results 
compatible with the ones obtained by using ICMMM2D within 
a relative error of 1%. The target precision of MMM2D in both 
cases has been set to 0.1%. A depletion layer due to entropic 
reasons is observed at the right interface, where no dielectric 
contrast is present. The opposite phenomenon is observed at 
the left interface, where the presence of induced surface charges 
introduces an effective attractive force which compensates the 
entropic depletion layer and generates a local density increase. 



4. Potential of Mean Force of a DNA segment in a Nanopore 

Experiments on DNA translocation through biological and 
synthetic pores have recently attracted a lot of attention. De- 
tailed descriptions of experiments are available e.g. in refs. [26, 
27, 28, 29]. It has been shown that the translocation rate of 
DNA through a pore depends strongly on the ionic strength of 
the buffer, hence indicating an electrostatic contribution to the 
translocation free energy barrier. Since a highly charged object 
like DNA is repelled from walls that are less polarizable than 
water, we investigate to which extent the DBFs influence the 
translocation free energy barrier. 

We modelled a synthetic nanopore using hard walls with a 
dielectric constant of 2, choosing the pore diameter and length 
to be 5 and 8 nm, respectively, as they fall within the range 
of experimental samples. On this length scale double stranded 
DNA (dsDNA) is stiff, as it has a persistence length of =* 50 nm 
at physiological conditions. Therefore it is justified to model 
it as a series of beads (500, in this work) constrained on the 
pore axis. The inter-bead distance is kept fixed, and each of 
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Figure 2: Free energy profiles for a DNA segment in the salt free case (solid 
lines, c = mmol/1) and with a salt concentration c = 100 mmol/1 (dashed 
lines). Both cases are reported with (squares) and without (circles) taking into 
account polarization effects with ICC*. Statistical error bars are smaller than 
the symbols. 



the beads bears a charge so that the linear charge density of 
0.17e/nm, characteristic of dsDNA, is reproduced. 

Ions and counterions are represented explicitly. They inter- 
act mutually and with the DNA beads through a WCA poten- 
tial having <x = 0.425 and <x = 1.215 nm, respectively, corre- 
sponding to a DNA diameter of 2 nm. The WCA interaction 
strength is set for both cases to 1 kgT. The electrostatic in- 
teractions are calculated in full 3D periodicity with the P3M 
algorithm. We employ a cubic simulation box with a 20 nm 
long edge. To mimic the case of low salt concentration, only 
monovalent counterions were introduced in the box so that the 
DNA molecule is neutralized. In an additional simulation, a 
finite salt concentration of 100 mmol/1 was added. All simula- 
tions were performed both with and without application of the 
ICC* algorithm to investigate the influence of DBFs. 

It is straightforward to calculate then the free energy barrier 
by computing the potential of mean force (PMF) acting on the 
center of mass of the model DNA along the pore axis. For this 
reaction coordinate the Fixman potential [30] is constant, and 
the PMF can be obtained by numerical integration of the mean 
force. The obtained PMFs are shown in Fig. 2. The free en- 
ergy barrier in the salt-free case is strikingly higher (increasing 
to approximately IQksT) when DBFs are taking into account 
using ICC*, in comparison to the case when DBFs are not con- 
sidered. On the contrary, at a salt concentration of 100 mmol/1, 
the barrier increase is less pronounced, and the curves obtained 
with or without taking into account DBFs show a comparable 
pattern with a barrier height of about 4ksT. The presence of a 
barrier for the higher salt concentration case, as well as in the 
case when no DBFs are considered, can be explained by the 
steric confinement of the counterion cloud in the nanopore. On 
the contrary, at low salt concentration the Coulomb interaction 
is not screened, and the effect of DBFs is maximized, leading 
to the observed higher free energy barrier. 



5. Conclusion 

We have presented the ICC* algorithm and have shown that 
the presence of dielectric boundary forces in coarse-grained 
simulations cannot be neglected under many conditions. Polar- 
ization effects due to the presence of dielectric mismatches at 
interfaces can lead to important effects, as it has been demon- 
strated for the case of a salt solution confined to a slit pore, 
as well as in the case of the translocation free energy barrier 
of a model DNA through a nanopore. Although for the DNA- 
nanopore system screening at finite salt concentration reduces 
the influence of the dielectric boundary forces, in the case of a 
slit pore the effects cannot be neglected even at relatively high 
(200 mmol/1) concentrations. 
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